﻿using System;

namespace MagistrAPI.MathModel
{
    public static class MathModelFunc
    {
        public static double H(double x, double y, double z)
        {
            return Rr(x, y, z) - R(Phiuu(x, y, z));
        }

        public static double Rr(double x, double y, double z)
        {
            return Math.Sqrt(x*x + y*y + z*z);
        }

        public static double R(double phiuu)
        {
            return MathModelConst.Rg*(1 - MathModelConst.Alphag*Math.Pow(Math.Sin(phiuu), 2));
        }

        public static double Phiuu(double x, double y, double z)
        {
            return Math.Asin(z/Rr(x, y, z));
        }

        public static double[,] Mrotv(double[,] matr, double e)
        {
            var s = 10e10;
            var k = 0;
            var v = new double[,] {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}};
            double[,] u = new double[3,3],
                      Jap,
                      c;
            double fi01, fi02, fi12;
            while (s > e)
            {
                s = matr[0, 1] + matr[0, 2] + matr[1, 2];
                if (Math.Abs(s) < e) break;
                if (Math.Abs(matr[0, 0] - matr[1, 1]) > e)
                    fi01 = 0.5*Math.Atan(-2*matr[0, 1]/(matr[0, 0] - matr[1, 1]));
                else fi01 = 0.5*Math.Atan(-2*matr[0, 1]/e);
                u[0, 0] = Math.Cos(fi01);
                u[0, 1] = Math.Sin(fi01);
                u[1, 0] = -Math.Sin(fi01);
                u[1, 1] = Math.Cos(fi01);
                u[2, 2] = 1;
                c = MatrToMatr(MatrT(u, 3, 3), matr, 3, 3, 3);
                Jap = MatrToMatr(c, u, 3, 3, 3);
                c = MatrToMatr(v, u, 3, 3, 3);
                v = c;
                if (Math.Abs(matr[0, 0] - matr[2, 2]) > e) fi02 = 0.5*Math.Atan(-2*Jap[0, 2]/(Jap[0, 0] - Jap[2, 2]));
                else fi02 = 0.5*Math.Atan(-2*matr[0, 2]/e);
                u = new double[3,3];
                u[0, 0] = Math.Cos(fi02);
                u[0, 2] = Math.Sin(fi02);
                u[1, 1] = 1;
                u[2, 0] = -Math.Sin(fi02);
                u[2, 2] = Math.Cos(fi02);
                c = MatrToMatr(MatrT(u, 3, 3), Jap, 3, 3, 3);
                Jap = MatrToMatr(c, u, 3, 3, 3);
                c = MatrToMatr(v, u, 3, 3, 3);
                v = c;
                if (Math.Abs(matr[1, 1] - matr[2, 2]) > e) fi12 = 0.5*Math.Atan(-2*Jap[1, 2]/(Jap[1, 1] - Jap[2, 2]));
                else fi12 = 0.5*Math.Atan(-2*matr[1, 2]/e);
                u = new double[3,3];
                u[0, 0] = 1;
                u[1, 1] = Math.Cos(fi12);
                u[1, 2] = Math.Sin(fi12);
                u[2, 1] = -Math.Sin(fi12);
                u[2, 2] = Math.Cos(fi12);
                c = MatrToMatr(MatrT(u, 3, 3), Jap, 3, 3, 3);
                Jap = MatrToMatr(c, u, 3, 3, 3);
                c = MatrToMatr(v, u, 3, 3, 3);
                v = c;
                matr = Jap;
                k++;
            }
            return v;
        }

        public static double[,] MatrToMatr(double[,] matrA, double[,] matrB, int nA, int mA, int mB)
        {
            var matrC = new double[nA,mB];
            double cur;
            for (int i = 0; i < nA; i++)
            {
                for (int k = 0; k < mB; k++)
                {
                    cur = 0;
                    for (int j = 0; j < mA; j++)
                    {
                        cur += matrA[i, j] * matrB[j, k];
                    }
                    matrC[i, k] = cur;
                }
            }
            return matrC;
        }

        public static double[,] MatrT(double[,] matr, int n, int m)
        {
            double cur;
            for (int i = 0; i < n - 1; i++)
            {
                for (int j = i + 1; j < m; j++)
                {
                    cur = matr[i, j];
                    matr[i, j] = matr[j, i];
                    matr[j, i] = cur;
                }
            }
            return matr;
        }

        public static double[,] LsbUy(double lX, double lY, double lZ, double mX, double mY, double mZ, double nX, double nY, double nZ)
        {
            var lmn = new[,] {{lX, lY, lZ}, {mX, mY, mZ}, {nX, nY, nZ}};
            return MathModelFunc.MatrToMatr(lmn, MathModelFunc.MatrT(MathModelMatrices.Jv, 3, 3), 3, 3, 3);
        }

        public static double[] Rxyz(double x, double y, double z, double vX, double vY, double vZ, double lX, double lY, double lZ,
                                    double mX, double mY, double mZ, double nX, double nY, double nZ)
        {
            double[] vbu = Vbu(x, y, z, vX, vY, vZ);
            double curRezult = MathModelConst.Sa*Q(Math.Sqrt(vbu[0]*vbu[0] + vbu[1]*vbu[1] + vbu[2]*vbu[2]), H(x, y, z));
            double rx = Cx(x, y, z, vX, vY, vZ, lX, lY, lZ, mX, mY, mZ, nX, nY, nZ)*curRezult;
            double ry = Cy(x, y, z, vX, vY, vZ, lX, lY, lZ, mX, mY, mZ, nX, nY, nZ)*curRezult;
            double rz = Cz(x, y, z, vX, vY, vZ, lX, lY, lZ, mX, mY, mZ, nX, nY, nZ)*curRezult;
            return new[] {rx, ry, rz};
        }

        public static double Q(double v, double rF)
        {
            return 0.5*Ph(rF)*v*v;
        }

        public static double Ph(double h)
        {
            return Math.Pow(Math.E, Lnp(h))/(MathModelConst.R0*Tm(h));
        }

        public static double Lnp(double h)
        {
            double g0 = MathModelConst.G0;
            double p0 = MathModelConst.P0;
            double r0 = MathModelConst.R0;
            double lNp11 = Math.Log(p0, Math.E) - g0/(r0*-0.00651122)*Math.Log(Tm(11000)/Tm(0), Math.E);
            double lNp25 = lNp11 - g0/(r0*Tm(25000))*(F(25000) - F(11000));
            double lNp46 = lNp25 - g0/(r0*0.00276098)*Math.Log(Tm(46000)/Tm(25000), Math.E);
            double lNp54 = lNp46 - g0/(r0*Tm(54000))*(F(54000) - F(46000));
            double lNp80 = lNp54 - g0/(r0*(-0.00349544))*Math.Log(Tm(80000)/Tm(54000), Math.E);
            double lNp95 = lNp80 - g0/(r0*Tm(95000))*(F(95000) - F(80000));
            double lNp110 = lNp95 - g0/(r0*(0.005))*Math.Log(Tm(110000)/Tm(95000), Math.E);
            double lNp120 = lNp110 - g0/(r0*Tm(120000))*(F(120000) - F(110000));
            double lNp150 = lNp120 - g0/(r0*(0.02346357))*Math.Log(Tm(150000)/Tm(120000), Math.E);
            double lNp160 = lNp150 - g0/(r0*Tm(160000))*(F(160000) - F(150000));
            if (h <= 11000) return Math.Log(p0, Math.E) - g0/(r0*(-0.00651122))*Math.Log(Tm(h)/Tm(0), Math.E);
            if (h > 11000 && h <= 25000) return lNp11 - g0/(r0*Tm(h))*(F(h) - F(11000));
            if (h > 25000 && h <= 46000) return lNp25 - g0/(r0*(0.00276098))*Math.Log(Tm(h)/Tm(25000), Math.E);
            if (h > 46000 && h <= 54000) return lNp46 - g0/(r0*Tm(h))*(F(h) - F(46000));
            if (h > 54000 && h <= 80000) return lNp54 - g0/(r0*(-0.00349544))*Math.Log(Tm(h)/Tm(54000), Math.E);
            if (h > 80000 && h <= 95000) return lNp80 - g0/(r0*Tm(h))*(F(h) - F(80000));
            if (h > 95000 && h <= 110000) return lNp95 - g0/(r0*(0.005))*Math.Log(Tm(h)/Tm(95000), Math.E);
            if (h > 110000 && h <= 120000) return lNp110 - g0/(r0*Tm(h))*(F(h) - F(110000));
            if (h > 120000 && h <= 150000) return lNp120 - g0/(r0*(0.02346357))*Math.Log(Tm(h)/Tm(120000), Math.E);
            if (h > 150000 && h <= 160000) return lNp150 - g0/(r0*Tm(h))*(F(h) - F(150000));
            if (h > 160000) return lNp160 - g0/(r0*(0.00308461))*Math.Log(Tm(h)/Tm(160000), Math.E);
            return 0;
        }

        public static double Tm(double h)
        {
            if (h <= 11000) return 288.15 - 0.00651122*F(h);
            if (h > 11000 && h <= 25000) return 216.66;
            if (h > 25000 && h <= 46000) return 216.66 + 0.00276098*(F(h) - F(25000));
            if (h > 46000 && h <= 54000) return 274;
            if (h > 54000 && h <= 80000) return 274 - 0.00349544*(F(h) - F(54000));
            if (h > 80000 && h <= 95000) return 185;
            if (h > 95000 && h <= 110000) return 185 + 0.005*(F(h) - F(95000));
            if (h > 110000 && h <= 120000) return 257.64 + 0.00801741*(F(h) - F(110000));
            if (h > 120000 && h <= 150000) return 335 + 0.02346357*(F(h) - F(120000));
            if (h > 150000 & h <= 160000) return 1010 + 0.01987408*(F(h) - F(150000));
            if (h > 160000) return 1199.4 + 0.00308461*(F(h) - F(160000));
            return 0;
        }

        public static double F(double h)
        {
            return h*(MathModelConst.Rz/(MathModelConst.Rz + h));
        }

        public static double[] Vbu(double x, double y, double z, double vX, double vY, double vZ)
        {
            double[] vect = Vect(0, 0, MathModelConst._Ω, x, y, z);
            return new[] {vX - vect[0], vY - vect[1], vZ - vect[2]};
        }

        public static double[] Vect(double x1, double y1, double z1, double x2, double y2, double z2)
        {
            return new[] {y1*z2 - y2*z1, z1*x2 - x1*z2, x1*y2 - x2*y1};
        }

        public static double Cx(double x, double y, double z, double vX, double vY, double vZ, double lX, double lY, double lZ,
                                double mX, double mY, double mZ, double nX, double nY, double nZ)
        {
            return MathModelConst.Cx0*Math.Cos(ααnb(x, y, z, vX, vY, vZ, lX, lY, lZ, mX, mY, mZ, nX, nY, nZ));
        }

        public static double ααnb(double x, double y, double z, double vX, double vY, double vZ, double lX, double lY, double lZ,
                                  double mX, double mY, double mZ, double nX, double nY, double nZ)
        {
            double[] vbsby = Vbsby(x, y, z, vX, vY, vZ, lX, lY, lZ, mX, mY, mZ, nX, nY, nZ);
            return Math.Acos(vbsby[0]/Math.Sqrt(vbsby[0]*vbsby[0] + vbsby[1]*vbsby[1] + vbsby[2]*vbsby[2]));
        }

        public static double[] Vbsby(double x, double y, double z, double vX, double vY, double vZ, double lX, double lY, double lZ,
                                     double mX, double mY, double mZ, double nX, double nY, double nZ)
        {
            double[] vbu = Vbu(x, y, z, vX, vY, vZ);
            double[,] vbsby = MathModelFunc.MatrToMatr(Luysb(lX, lY, lZ, mX, mY, mZ, nX, nY, nZ),
                                                       new[,] {{vbu[0]}, {vbu[1]}, {vbu[2]}}, 3, 3, 1);
            return new[] {vbsby[0, 0], vbsby[1, 0], vbsby[2, 0]};
        }

        public static double[,] Luysb(double lX, double lY, double lZ, double mX, double mY, double mZ, double nX, double nY, double nZ)
        {
            var lmn = new[,] {{lX, mX, nX}, {lY, mY, nY}, {lZ, mZ, nZ}};
            return MathModelFunc.MatrToMatr(MathModelMatrices.Jv, lmn, 3, 3, 3);
        }

        public static double Cyn(double x, double y, double z, double vX, double vY, double vZ, double lX, double lY, double lZ,
                                 double mX, double mY, double mZ, double nX, double nY, double nZ)
        {
            return MathModelConst.Cyn0*Math.Sin(ααnb(x, y, z, vX, vY, vZ, lX, lY, lZ, mX, mY, mZ, nX, nY, nZ));
        }

        public static double Cy(double x, double y, double z, double vX, double vY, double vZ, double lX, double lY, double lZ,
                                double mX, double mY, double mZ, double nX, double nY, double nZ)
        {
            double[] vbsby = Vbsby(x, y, z, vX, vY, vZ, lX, lY, lZ, mX, mY, mZ, nX, nY, nZ);

            return Cyn(x, y, z, vX, vY, vZ, lX, lY, lZ, mX, mY, mZ, nX, nY, nZ)*
                   ((-vbsby[1])/Math.Sqrt(vbsby[1]*vbsby[1] + vbsby[2]*vbsby[2]));
        }

        public static double Cz(double x, double y, double z, double vX, double vY, double vZ, double lX, double lY, double lZ,
                                double mX, double mY, double mZ, double nX, double nY, double nZ)
        {
            double[] vbsby = Vbsby(x, y, z, vX, vY, vZ, lX, lY, lZ, mX, mY, mZ, nX, nY, nZ);

            return -Cyn(x, y, z, vX, vY, vZ, lX, lY, lZ, mX, mY, mZ, nX, nY, nZ)*
                   (vbsby[2]/Math.Sqrt(vbsby[1]*vbsby[1] + vbsby[2]*vbsby[2]));
        }

        public static double[] Guu(double x, double y, double z)
        {
            var lLfu = new[,]
                           {
                               {Math.Cos(φuu(x, y, z)), 0, Math.Sin(φuu(x, y, z))},
                               {0, 1, 0},
                               {-Math.Sin(φuu(x, y, z)), 0, Math.Cos(φuu(x, y, z))}
                           };
            var lLgamu = new[,]
                             {
                                 {Math.Cos(λλu(x, y, z)), Math.Sin(λλu(x, y, z)), 0},
                                 {-Math.Sin(λλu(x, y, z)), Math.Cos(λλu(x, y, z)), 0},
                                 {0, 0, 1}
                             };
            var gg = new double[3,1] {{Gr(x, y, z)}, {Gφ(x, y, z)}, {Gλ(x, y, z)}};
            var guu =
                MathModelFunc.MatrToMatr(
                    MathModelFunc.MatrToMatr(MathModelFunc.MatrT(lLgamu, 3, 3), MathModelFunc.MatrT(lLfu, 3, 3), 3, 3,
                                             3), gg, 3, 3, 1);
            return new[] {guu[0, 0], guu[1, 0], guu[2, 0]};
        }

        public static double φuu(double x, double y, double z)
        {
            return Math.Asin(z/Math.Sqrt(x*x + y*y + z*z));
        }

        public static double Arcsin(double sin, double cos)
        {
            double x = Math.Asin(sin);
            double z;
            if (cos < 0)
            {
                if (sin < 0)
                {
                    z = -(x + Math.PI);
                }
                else
                {
                    z = -(x - Math.PI);
                }
            }
            else
            {
                z = x;
            }
            return z;
        }

        public static double λλu(double x, double y, double z)
        {
            return Arcsin(y/Math.Sqrt(x*x + y*y), x/Math.Sqrt(x*x + y*y));
        }

        public static double AbsV(double x, double y, double z)
        {
            return Math.Sqrt(x*x + y*y + z*z);
        }

        public static double Gr(double x, double y, double z)
        {
            double pi2 = -1.756 * Math.Pow(10, 25);
            double pi4 = 1.548 * Math.Pow(10, 36);
            double ff = φuu(x, y, z);
            double gr = 1 / Math.Pow(AbsV(x, y, z), 2) *
                        (-MathModelConst.K + 1.5 * pi2 / Math.Pow(AbsV(x, y, z), 2) * (1 - 3 * Math.Pow(ff, 2)) -
                         5 / 8.0 * pi4 / Math.Pow(AbsV(x, y, z), 4) *
                         (35 * Math.Pow(Math.Sin(ff), 4) - 30 * Math.Pow(Math.Sin(ff), 2) + 3));
            return gr;
        }

        public static double Gφ(double x, double y, double z)
        {
            double pi2 = -1.756 * Math.Pow(10, 25);
            double pi4 = 1.548 * Math.Pow(10, 36);
            double ff = φuu(x, y, z);
            double gf = 1 / Math.Pow(AbsV(x, y, z), 4) *
                           (3 * pi2 + 2.5 * pi4 / Math.Pow(AbsV(x, y, z), 2) * (7 * Math.Pow(Math.Sin(ff), 2) - 3)) * Math.Sin(ff) *
                           Math.Cos(ff);
            return gf;
        }

        public static double Gλ(double x, double y, double z)
        {
            return 0;
        }

        public static double Mx(double x, double y, double z, double vX, double vY, double vZ, double lX, double lY, double lZ,
                                double mX, double mY, double mZ, double nX, double nY, double nZ)
        {
            double[] vbu = Vbu(x, y, z, vX, vY, vZ);
            double qY = Q(AbsV(vbu[0], vbu[1], vbu[2]), H(x, y, z));
            return (Mxx(x, y, z, vX, vY, vZ, lX, lY, lZ, mX, mY, mZ, nX, nY, nZ) +
                    MathModelConst.DeltaZ*Cy(x, y, z, vX, vY, vZ, lX, lY, lZ, mX, mY, mZ, nX, nY, nZ) -
                    MathModelConst.DeltaY*Cz(x, y, z, vX, vY, vZ, lX, lY, lZ, mX, mY, mZ, nX, nY, nZ))*MathModelConst.Sa*
                   qY;
        }

        public static double My(double x, double y, double z, double vX, double vY, double vZ, double lX, double lY, double lZ,
                                double mX, double mY, double mZ, double nX, double nY, double nZ)
        {
            double[] vbu = Vbu(x, y, z, vX, vY, vZ);
            double qY = Q(AbsV(vbu[0], vbu[1], vbu[2]), H(x, y, z));
            return (Myy(x, y, z, vX, vY, vZ, lX, lY, lZ, mX, mY, mZ, nX, nY, nZ) +
                    MathModelConst.DeltaX*Cz(x, y, z, vX, vY, vZ, lX, lY, lZ, mX, mY, mZ, nX, nY, nZ) -
                    MathModelConst.DeltaZ*Cx(x, y, z, vX, vY, vZ, lX, lY, lZ, mX, mY, mZ, nX, nY, nZ))*MathModelConst.Sa*
                   qY*MathModelConst.La;
        }

        public static double Mz(double x, double y, double z, double vX, double vY, double vZ, double lX, double lY, double lZ,
                                double mX, double mY, double mZ, double nX, double nY, double nZ)
        {
            double[] vbu = Vbu(x, y, z, vX, vY, vZ);
            double qY = Q(AbsV(vbu[0], vbu[1], vbu[2]), H(x, y, z));
            return (Mzz(x, y, z, vX, vY, vZ, lX, lY, lZ, mX, mY, mZ, nX, nY, nZ) +
                    MathModelConst.DeltaY*Cx(x, y, z, vX, vY, vZ, lX, lY, lZ, mX, mY, mZ, nX, nY, nZ) -
                    MathModelConst.DeltaX*Cy(x, y, z, vX, vY, vZ, lX, lY, lZ, mX, mY, mZ, nX, nY, nZ))*MathModelConst.Sa*
                   qY*MathModelConst.La;

        }

        public static double Mzn(double x, double y, double z, double vX, double vY, double vZ, double lX, double lY, double lZ,
                                 double mX, double mY, double mZ, double nX, double nY, double nZ)
        {
            return -MathModelConst.Delta*MathModelConst.Cyn0*
                   Math.Sin(ααnb(x, y, z, vX, vY, vZ, lX, lY, lZ, mX, mY, mZ, nX, nY, nZ));
        }

        private static double Mxx(double x, double y, double z, double vX, double vY, double vZ, double lX, double lY, double lZ,
                                 double mX, double mY, double mZ, double nX, double nY, double nZ)
        {
            return 0;
        }

        private static double Myy(double x, double y, double z, double vX, double vY, double vZ, double lX, double lY, double lZ,
                                 double mX, double mY, double mZ, double nX, double nY, double nZ)
        {
            double[] vbsby = Vbsby(x, y, z, vX, vY, vZ, lX, lY, lZ, mX, mY, mZ, nX, nY, nZ);
            return Mzn(x, y, z, vX, vY, vZ, lX, lY, lZ, mX, mY, mZ, nX, nY, nZ) *
                   (vbsby[2] / Math.Sqrt(vbsby[1] * vbsby[1] + vbsby[2] * vbsby[2]));
        }

        private static double Mzz(double x, double y, double z, double vX, double vY, double vZ, double lX, double lY, double lZ,
                                 double mX, double mY, double mZ, double nX, double nY, double nZ)
        {
            double[] vbsby = Vbsby(x, y, z, vX, vY, vZ, lX, lY, lZ, mX, mY, mZ, nX, nY, nZ);
            return Mzn(x, y, z, vX, vY, vZ, lX, lY, lZ, mX, mY, mZ, nX, nY, nZ)*
                   (-vbsby[1]/Math.Sqrt(vbsby[1]*vbsby[1] + vbsby[2]*vbsby[2]));
        }

        public static double[,] LsbTr(double fi, double al, double gam)
        {
            double[,] curRezult = MatrToMatr(MathModelMatrices.Lφn(fi), MathModelMatrices.Lαn(al), 3, 3, 3);
            return MatrT(MatrToMatr(curRezult, MathModelMatrices.Lγn(gam), 3, 3, 3), 3, 3);
        }

        public static double[,] Ltru(double x, double y, double z, double vX, double vY, double vZ)
        {
            return new[,]
                       {
                           {Ikx(vX, vY, vZ), Jkx(x, y, z, vX, vY, vZ), Kkx(x, y, z, vX, vY, vZ)},
                           {Iky(vX, vY, vZ), Jky(x, y, z, vX, vY, vZ), Kky(x, y, z, vX, vY, vZ)},
                           {Ikz(vX, vY, vZ), Jkz(x, y, z, vX, vY, vZ), Kkz(x, y, z, vX, vY, vZ)}
                       };
        }

        private static double Ikx(double vX, double vY, double vZ)
        {
            return vX/AbsV(vX, vY, vZ);
        }

        private static double Iky(double vX, double vY, double vZ)
        {
            return vY / AbsV(vX, vY, vZ);
        }

        private static double Ikz(double vX, double vY, double vZ)
        {
            return vZ / AbsV(vX, vY, vZ);
        }

        private static double Jkx(double x, double y, double z, double vX, double vY, double vZ)
        {
            double[] vect = Vect(Kkx(x, y, z, vX, vY, vZ), Kky(x, y, z, vX, vY, vZ), Kkz(x, y, z, vX, vY, vZ),
                                 Ikx(vX, vY, vZ), Iky(vX, vY, vZ), Ikz(vX, vY, vZ));
            return vect[0];
        }

        private static double Jky(double x, double y, double z, double vX, double vY, double vZ)
        {
            double[] vect = Vect(Kkx(x, y, z, vX, vY, vZ), Kky(x, y, z, vX, vY, vZ), Kkz(x, y, z, vX, vY, vZ),
                                 Ikx(vX, vY, vZ), Iky(vX, vY, vZ), Ikz(vX, vY, vZ));
            return vect[1];
        }

        private static double Jkz(double x, double y, double z, double vX, double vY, double vZ)
        {
            double[] vect = Vect(Kkx(x, y, z, vX, vY, vZ), Kky(x, y, z, vX, vY, vZ), Kkz(x, y, z, vX, vY, vZ),
                                 Ikx(vX, vY, vZ), Iky(vX, vY, vZ), Ikz(vX, vY, vZ));
            return vect[2];
        }

        private static double Kkx(double x, double y, double z, double vX, double vY, double vZ)
        {
            double[] guu = Guu(x, y, z);
            double[] vect = Vect(guu[0], guu[1], guu[2], vX, vY, vZ);
            return vect[0]/AbsV(vect[0], vect[1], vect[2]);
        }

        private static double Kky(double x, double y, double z, double vX, double vY, double vZ)
        {
            double[] guu = Guu(x, y, z);
            double[] vect = Vect(guu[0], guu[1], guu[2], vX, vY, vZ);
            return vect[1] / AbsV(vect[0], vect[1], vect[2]);
        }

        private static double Kkz(double x, double y, double z, double vX, double vY, double vZ)
        {
            double[] guu = Guu(x, y, z);
            double[] vect = Vect(guu[0], guu[1], guu[2], vX, vY, vZ);
            return vect[2] / AbsV(vect[0], vect[1], vect[2]);
        }

        public static double[,] LugL0(double x, double y, double z, double vX, double vY, double vZ, double fi, double al, double gam)
        {
            double[] vbu = Vbu(x, y, z, vX, vY, vZ);
            double[,] lsbu = MatrToMatr(Ltru(x, y, z, vbu[0], vbu[1], vbu[2]), LsbTr(fi, al, gam), 3, 3, 3);
            return MatrToMatr(MatrT(MathModelMatrices.Jv, 3, 3), MatrT(lsbu, 3, 3), 3, 3, 3);
        }
    }
}